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Using in situ low-energy electron microscopy and density functional theory, we studied the growth 
structure and work function of bilayer graphene on Pd(lll). Low-energy electron diffraction analysis 
established that the two graphene layers have multiple rotational orientations relative to each other 
and the substrate plane. We observed heterogeneous nucleation and simultaneous growth of multiple, 
faceted layers prior to the completion of second layer. We propose that the facetted shapes are 
due to the zigzag-terminated edges bounding graphene layers growing under the larger overlying 
layers. We also found that the work functions of bilayer graphene domains are higher than those of 
monolayer graphene, and depend sensitively on the orientations of both layers with respect to the 
substrate. Based on first-principles simulations, we attribute this behavior to oppositely oriented 
electrostatic dipoles at the graphene/Pd and graphene/graphene interfaces, whose strengths depend 
on the orientations of the two graphene layers. 



I. INTRODUCTION 

Few-layer graphene^ is attractive for applications in 
optoelectronics as transparent conductors -where precise 
control over layer thickness is not necessary- owing to 
its low sheet resistance and high transparency^ and 
in field-effect transistors due to the opening of an elec- 
tronic bandgap that, at least in bilayer graphene, can be 
tuned with an electric field. ^ Given that device charac- 
teristics depend on the work function of graphene and 
on how it varies at contacts,^ it is of fundamental and 
practical importance to understand the nature of the 
contact (Ohmic or Schottky) at metal-graphene inter- 
faces. In the case of monolayer graphene, recent studies 
indicate that the crystallographic orientation of graphene 
domains with respect to the metal can alter their elec- 
tronic properties. For few-layer graphene, the role of 
the specific metal surface and of the in-plane orienta- 
tions of those layers on the electrical characteristics of 
the graphene is not well known. We choose Pd(lll), one 
of the commonly used contact materials recently shown 
to exhibit low contact resistance,^ as a model to inves- 
tigate the influence of thickness and orientation on the 
work functions of graphene layers on metals. 

Here, we present results of in situ low-energy elec- 
tron microscopy (LEEM) and density functional theory 
(DFT) investigations of the growth structure and work 
function of bilayer graphene on Pd(lll). Selected area 
low-energy electron diffraction (LEED) patterns indicate 
that Bernal stacking is typically not observed in the as- 
grown graphene layers. From the electron reflectivity 
data, we have extracted work functions of graphene do- 



mains as a function of the in-plane orientations of the 
top and bottom layers. For monolayer graphene, we have 
previously shown that the work functions can vary by up 
to 0.15 eV depending on the orientation of the domains 
with respect to the substrate^^ Different registries be- 
tween the graphene monolayer and the Pd(lll) substrate 
lead to strong spatial variations of the charge transfer 
between the graphene and substrate. The resulting net 
surface dipole moment varies with the in-plane orienta- 
tion of the monolayer. In case of bilayer graphene, we 
find that a smaller net electric dipole develops between 
the topmost graphene layer and the rest of the system 
(i.e., the monolayer-covered substrate). The direction 
of this secondary dipole moment is opposite to that of 
the dipole developed at the Pd/graphene interface, and 
leads to an increase in the work function as compared 
to monolayer-covered Pd. Furthermore, the magnitude 
of the secondary dipole changes with the orientation of 
the second layer. This observation suggests that the first 
graphene layer does not fully passivate the substrate but 
allows charge transfer into the second layer, thus affecting 
the work function of the supported bilayer graphene. 



II. EXPERIMENTAL METHODS AND 
RESULTS 

A. Methods 

All of our experiments are carried out using a carbon- 
saturated Pd(lll) single crystal in an ultra-high vac- 
uum (UHV, base pressure < 1.0 x 10"^° Torr) LEEM 
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system.^'^ Sample preparation and other experimental 
details are presented in Ref. ?. Graphene layers of desired 
thickness are obtained by cooling from ~ 900 °C to lower 
temperatures, during which C segregates from the bulk 
to the surface. ^^^^ For example, monolayer graphene is 
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obtained when the sample is cooled from ~ 900 to ~ 790 
°C at the rate of ~ 1 K/s. Upon completion of the first 
layer, further cooling to ~ 700 °C or lower yields multi- 
layer graphene. Bright-field LEEM images are acquired 
as a function of time at the annealing temperature. The 
direct observation of graphene formation helps identify 
the layer thickness and is also useful in understanding 
the nucleation and growth kinetics. After the growth of 
a desired number of layers, the sample is quench-cooled to 
room temperature. LEEM images are then acquired as a 
function of electron energy E at 0.1 eV intervals between 
E = -5 and 30 eV and 0.01 eV intervals between E = 
0.0 and 2.5 eV. From the LEEM image intensity I ys. E 
data, we determine the graphene layer thickness, the elec- 
tron injection threshold energy ^, and the work function 
Selected area LEED patterns^ are used to iden- 
tify the orientation 6 of graphene layers with respect to 
the substrate. We define as the angle between Pd[110] 
and graphene [11 20] directions with the positive (nega- 
tive) sign denoting in-plane, clockwise (counterclockwise) 
rotation. The measurement uncertainties in values are 
±1°. 



B. LEEM results 




FIG. 1: Representative LEEM images of Pd(lll) acquired 
after cooling from 790 °C to 660 °C during the growth of 
multilayer graphene at times t = (a) 17 s, (b) 39 s, (c) 83 s, 
and (d) 205 s. (t = corresponds to the time at which the 
sample temperature is 660 °C.) Field of view = 9.3 /xm and 
electron energy E = 4.3 eV. 



Figure [T] shows a series of bright-field LEEM images 
acquired from a Pd(lll) sample during the growth of 
multilayer graphene at 660 °C. The sample is cooled from 
790 °C when the surface was already covered with 1 ML 



graphene (Fig. [T^). The alternating brighter and darker 
grey features in Figs. [TjD-d, as we show below, are multi- 
layer stacks of graphene. (The black contrast region ob- 
served along the upper right corner of the LEEM images 
is due to a three-dimensional Pd mound on the surface.) 
We observe three interesting phenomena: (i) growth oc- 
curs at selected locations on the sample, suggestive of 
heterogeneous nucleation, presumably at surface defects; 
(ii) multiple layers nucleate and grow simultaneously, i.e., 
third and higher layers appear before the completion of 
the second layer. In this measurement sequence, within 
205 s of cooling to 660 °C, most of the nucleated sites 
are covered with ten or more layers of graphene at the 
nucleation sites (see Fig.[l]i); and (iii) the multilayer do- 
mains are faceted with regular triangular and/or trun- 
cated hexagonal shapes, indicative of strong anisotropy 
in step energies. 

We explain this observation as follows. Graphene's 
honeycomb lattice gives rise to two types of simple edge 
structures, armchair and zigzag. The angle between any 
two armchair or any two zigzag edges is 60° or 120°, 
and that between armchair and zigzag edges is 30° or 
90°. Therefore, if the edges of the equilateral triangu- 
lar and/or truncated hexagonal shapes are simple, they 
must all be either zigzag or armchair in structure. By 
comparing the directions of domain edges observed in 
the LEEM images with the LEED spot orientations, we 
determined that the domains are bounded by graphene 
sheets with zigzag edges. We suggest that asymmetries 
in the geometry of zigzag edges relative to an adjacent 
graphene layer may account for the anisotropic shapes 
of multilayer domains. The asymmetry is most easily 
illustrated for the case of two Bernal-stacked layers, as 
shown in Fig.[2j All of the armchair edges are equivalent. 
However, any two zigzag edges oriented at 120° to each 
other are not equivalent, analogous to the A- and B- type 
steps on (lll)-oriented face-centered cubic crystals. As 
we describe in the next section, however, none of the ob- 
served graphene bilayer domains exhibit Bernal stacking 
(see Fig. [3]). Yet, it is possible that the energies of both 
armchair and zigzag edges can vary with their orienta- 
tion relative to the adjacent layer. This may explain the 
observation of multiple types of anisotropic shapes. 

In the following sections, we report the thickness, 
crystallographic orientation, and work functions of the 
graphene layers. Fig. [3^ is a representative LEEM im- 
age of multilayer graphene grown by cooling the sample 
from 900 °C to 709 °C and held for 1200 s. Fig. [3]d is 
a typical plot of / vs. E data collected from the regions 
(color coded for clarity) highlighted in Fig. [3^1. Note the 
oscillations in / at £^ values between ~ 1 and ^ 8 eV. 
These oscillations are a direct consequence of interference 
between electrons refiecting from the surface and those 
from interlayer and/or layer-substrate interfaces The 
number of oscillations increases with increasing number 
of layers. For graphene on Pd(lll), we find that n-layer 
graphene, irrespective of the interlayer stacking, exhibits 
(n-1) oscillations, i.e., 1, 2, 3, and 4 ML show 0, 1, 2, and 
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FIG. 2: Top- view schematic of Bernal- stacked bilayer 
graphene lattice. The top and bottom panels show the upper 
graphene domain bounded by armchair and zigzag edges, re- 
spectively. For clarity, the lower graphene layer is wider and 
shows the carbon network, while the upper layer shows only 
the edge atoms and their first neighbors. 



3 oscillations, respectively (Fig.[3]3). This result is similar 
to the / vs. E data for multilayer graphene on SiC(OOOl) 
after accounting for the so-called "buffer" layer J^iii 



C. LEED results 

Next, we have determined the crystallographic orienta- 
tions of the individual graphene layers using selected area 
LEED in combination with / vs. E data. Figs. [3t-f are 
representative LEED patterns acquired from the regions, 
highlighted in Fig. [3^, of 1, 2, 3, and 4 ML graphene 
analyzed above. We first identify the diffraction spots 
corresponding to Pd(lll). Of all the LEED patterns in 
Figs. [3t-f, only Fig. [3t shows six-fold symmetric spots 
corresponding to the Pd(lll)-1 x 1 lattice. The other set 
of spots in Fig. [3t are due to the monolayer graphene- 
1x1 lattice, which is rotated with respect to the substrate 
by ^ = —14°. This data serves as the reference for the 
identification of LEED spots corresponding to 2, 3, and 
4 ML graphene in Figs.[3]i-f. Since the LEED spot inten- 
sity decreases with increasing distance into the surface, 
we can determine the position of the graphene layer with 
respect to the surface. 

For example, in Fig. [3]i, we observe two sets of six- 
fold symmetric spots due to bilayer graphene, with the 
stronger and weaker intensity spots oriented along 9 = 
-14° and = +14°, respectively. That is, the top 




FIG. 3: (a) Typical LEEM image of multilayer graphene on 
Pd(lll). Field of view = 9.3 /xm, E = 3.25 eV. (b) Image 
intensity / vs. E for multilayer graphene. (c-f) Selected area 
LEED patterns of 1, 2, 3, and 4 ML graphene-covered regions 
highlighted by blue, green, purple, and orange squares in (a), 
respectively. Red arrow in (c) indicates the Pd[110]. The blue 
arrow indicates the graphene [11 20] for the topmost-graphene 
sheet. The green and purple arrows indicate graphene [11 20] 
for the buried graphene sheets that give the second strongest 
diffraction from 2 and 3 ML, respectively. 



graphene layer is oriented along Oi = —14° while the 
layer closer to the substrate is along 62 = +14°. The 
LEED patterns from 3 and 4 ML graphene, in Figs. [3^ 
and[3f, respectively also show two sets of six-fold symmet- 
ric spots with the strongest spots along 6 = —14° and the 
weaker spots along = +8.9°. We identify = —14° as 
the topmost layer, in both Figs. 3e and 3f. However, we 
could not determine which of the underlying graphene 
layers gives rise to the = +8.9° spots. In all of the 
LEED patterns, we find that one of the layers in multi- 
layer graphene is oriented at the same angle ( = —14° 
in the example of Fig. 3) as in the 1 ML graphene that 
surrounds the multilayer stack. This observation sug- 
gests that orientation of the top layer graphene is unaf- 
fected during the growth of subsequent layers, consistent 
with the expected growth of multiple layers below the 
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FIG. 4: Experimental work function values for bilayer 
graphene on Pd(lll) for different orientations Oi and O2 of 
the two graphene layers. 



SiC (where the ^ of bilayer graphene is also higher than 
that of a monolayer ),^^ but opposite to those for free- 
standing graphene transferred onto Si02,^^ where ^ is 
found to decrease with increasing layer thickness. We 
note, however, that the effect of domain orientation on 
the ^ values of bilayer graphene on Pd(lll) is smaller 
(up to 0.07 eV) compared to that observed (up to 0.15 
eV) in <l> for monolayer graphene domains (see Fig. 4). 
In order to confirm these results, we determined ^ values 
using a different procedure, ^^^^^ and obtained consistent 
results. An important factor that affects the work func- 
tion of bilayer graphene is the in-plane orientation of the 
two layers, which we next discuss from the perspective of 
electronic structure calculations. 



III. DFT RESULTS AND DISCUSSION 



layer that grows first. From LEED patterns acquired 
from over 10 different regions of the sample, we identi- 
fied at least 12 different in-plane orientations in bilayer 
graphene. In all cases the two layers are not rotation- 
ally aligned, as required for Bernal stacking, similar to 
previous reports.^ Based upon our data, we suggest that 
Bernal stacking will be rare in graphene layers grown on 
Pd(lll). 



D. Work function of graphene bilayers on Pd 

We now focus on the relationship between work func- 
tion ^ and orientation of bilayer graphene. Following 
the procedure outhned in Ref. 7, we first measure the 
electron injection threshold energy 0, at which the im- 
age intensity / drops to 90% of the intensity value at E 
= eV. Then, $ = -\- cj)^ where ^fu is the work 
function of the electron gun filament. From the / vs. 
E data obtained from clean Pd(lll) along with known 
values of $ for Pd (5.3—5.6 eV),^ we estimate ^fu to 
be between 3.1 and 3.4 eV. Given the large uncertain- 
ties in the knowledge of ^ju (up to 0.3 eV), we do not 
emphasize determining absolute values of ^. However, (j) 
values are precise to within 0.02 eV and, therefore, can 
be used to compare how $ changes with orientation and 
layer thickness. Recently, we reported that ^ of mono- 
layer graphene on Pd(lll) varies with in-plane rotation 
due to spatial variations in charge transfer at graphene 
- Pd interface.^ Using the same approach and assuming 
that ^fii = 3.4 eV, we have determined Phi for bilayer 
graphene domains as a function of Oi and O2. 

Fig. 4 shows the work function values ^ for bilayer 
graphene plotted at different values of 61 and 62. For 
all domain orientations, the work function of bilayer 
graphene is higher than that for monolayer graphene. 
Our results for ^ of graphene bilayers on Pd(lll) are 
similar to recent reports on epitaxial graphene grown on 



In order to gain insights into how the second graphene 
layer and the orientation of both layers affect the work 
function, we have performed DFT calculations on sys- 
tems with one and two graphene layers placed on one side 
of a Pd(lll) substrate. We have used the Siesta code^^ 
with double-zeta (DZ) basis set within the local density 
approximation (LDA) with the Ceperley- Alder exchange- 
correlation functional. A uniform grid in real space with 
a mesh cutoff of 150 Ry was used. The Brillouin zone 
has been sampled only at the F point since the unit cells 
were large, with periodicity vectors of 19.72A. The resid- 
ual forces on any atom are smaller than 0.04 eV/A at the 
end of structural relaxations, and the total energy has 
been converged to within 10 ~^ eV for all electronic prop- 
erty calculations. We have used 8x8 surface super cells 
for the construction of 19.9° and 30°-rotated graphene 
domains (126 C atoms and 128 C atoms, respectively), 
and 46 Pd atoms per layer with four Pd layers in the 
substrate. For the 30° orientation, we have adopted the 
epitaxial structure proposed by Giovanetti et a/.^^ The 
choice of the two particular orientations, 19.9° and 30°, 
is determined both by the necessity to keep the computa- 
tions tractable and by the fact that these two orientations 
can be placed simultaneously on the Pd(lll) substrate 
with virtually no strain in either of them. Thus, any 
computed variations in $ are due only to their orienta- 
tions and/or stacking order, but not to the (insignificant) 
strain in the graphene layers. 

For monolayer graphene, we compute the charge trans- 
fer and the associated dipole moment in a manner simi- 
lar to previous reports, but strictly adopt the standard 
direction convention for electrostatic dipole (i.e., dipole 
vector pointing from the negative charge to the positive 
charge). For single- and bi-layer graphene, the planar av- 
eraged electron density transferred with respect to all the 
individual components (substrate, first layer, and second 
layer -when present), is defined by 

Apa{z) = PGr2/Gri/Pd(^) " PPd{z) - pGn (^) - PGr2(^)- 

(1) 
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TABLE I: Work function results from DFT calculations for single- 


and bi- layer graphene, along 


with the net dipole Pa and 


the secondary dipole Pb 


defined in text. 


The work function of pure 


Pd substrate is also given as 


a reference. The electronic 


transfer to the graphene 


layer(s) is given in the last column, where 1 


layer denotes the one closest to the Pd substrate. 
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FIG. 5: Work function values from dipole moment calcula- 
tions, $re/ — eP/eo, vs. values determined directly from DFT 
calculations, (a) Bare Pd surface is taken as the reference 
for single- and bi-layer epitaxial systems, (b) Pd substrate 
covered with monolayer graphene is used as reference for the 
bi-layer calculations. 



This electron transfer creates a net dipole moment when 
integrated from the middle of the substrate to the middle 
of the vacuum in the supercell, Pa = J zApa{z)dz^ dipole 
which should be responsible for the change in ^ with 
respect to the bare Pd surface 



^Gr2/Gri/Pd = ^Pd - ePa/cQ. 



(2) 



Indeed, computing the work function directly from the 
electrostatic potential output of the DFT simulation 
(= ^vacuum " ^^Fermi) and from the dipole moment Pa 
[via Eq. (2)] shows that the electrostatic dipole largely 
accounts for the work function of single- and bi-layer 
graphene on Pd (Fig. 5a). For both single- and bi-layer 
graphene, the values of the dipole moment Pa are pos- 
itive (i.e., pointing away from the substrate, see Table 
1) and of same order of magnitude as that reported for 
other graphene-metal systems. 

To account specifically for the influence of the second 
graphene layer on <l>, we resort to another way of com- 



puting the dipole moment, one in which the reference 
consists of the topmost (second) graphene layer and the 
palladium substrate covered with monolayer graphene, 

Apb{z) = pGr2/Gri/Pd(^) " PGri/Pd(^) " PGr^iz)- (3) 

The charge transfer Api){z) and the dipole moment P5 = 
J zApij{z)dz developed with respect to the new reference 
are, therefore, much smaller when compared to the trans- 
fer to the first layer, as shown in Fig. 6 and listed in Ta- 
ble 1. Although smaller, the secondary dipole moment 
P5 does account for the work function changes between 
the single- and bi- layer graphene on Pd, 



Gr2/Gri/I 



i/Pd = *^Gri/Pd - ePb/eo, 



(4) 



as seen in Fig. 4b. The secondary dipole moment P5 
(Table 1) is always negative (i.e., it points into the sub- 
strate), which explains our observation that the work 
functions of bilayer systems are greater than those of 
monolayer systems. We note that although the varia- 
tions in experimentally measured values are in qualita- 
tive agreement with the DFT variations, they are smaller 
than those obtained using DFT. Possible reasons for the 
smaller variations in experimental values are (a) differ- 
ent orientations used in calculations compared to exper- 
iments (our initial attempts to use the experimental ori- 
entations in DFT simulations resulted in supercells that 
were either too large or too strained), and (b) the pres- 
ence of elastic strain, carbon adatoms at the graphene-Pd 
interface, or other impurities in graphene layers. 

Finally, we have also estimated the electron transfer 
to each graphene layer by integrating the density using 
as integration limits the midpoints between atomic lay- 
ers. For computing the electron transfer to the second 
graphene layer (topmost), we integrated from the mid- 
point between the first and second graphene layers to 
the middle of the vacuum. The results, listed in the last 
column of Table 1, show that the first graphene layer 
is p-doped (loss of electrons), while the second layer is 
n-doped. As seen in Table 1, the magnitude of the elec- 
tron transfer to the second layer is relatively small and 
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FIG. 6: Plane- averaged electron density transferred (number of electrons per unit volume) and surface dipole densities for (a) 
30.0° /Pd and (b) 19.9°/30.0°/Pd. The vertical gray (blue) lines indicate the locations of the graphene (substrate) layers. The 
origin of the z coordinate is taken in the middle of the corresponding reference substrates, i.e., (a) the midpoint of the bare Pd 
slab, and (b) the midpoint of the segment between inner (first) graphene layer and bare Pd layer on the opposite side of the 
slab. The arrows in the insets show schematically (a) the surface dipole along the surface normal for monolayer graphene, and 
(b) the secondary dipole pointing opposite to the surface normal. 



can be attributed to charge screening by the first layer. 
The signs of the estimated electronic transfer to each 
graphene layer are certainly consistent with directions of 
the calculated dipoles Pa and P5, i.e., a loss of electrons 
on the first graphene layer leads to a dipole Pa parallel 
to the surface normal, while a gain of electrons on the 
second graphene layers leads to a secondary dipole P5 
antiparallel to the surface normal. 

IV. CONCLUSIONS 

In conclusion, we used in situ LEEM and DFT to inves- 
tigate the growth structure and work function of multi- 
layer graphene on Pd(lll). None of the bilayer domains 
we examined had the rotational alignment required for 
Bernal stacking. Instead, there are multiple orientations 
of the layers relative to the substrate and each other. The 
orientation-dependent variations in work functions of bi- 



layer graphene are relatively smaller than those observed 
in monolayer graphene. We explain these results using 
DFT calculations, which reveal that the charge transfer 
from Pd depends sensitively on the orientations of both 
graphene layers. Our results indicate that one cannot 
infer how multilayers of graphene interact with a metal 
substrate (for example, in a contact) by examining only 
the first layer. Therefore, for the fabrication of graphene- 
based transistors and/or transparent conductors, knowl- 
edge of the graphene layer stacking and its influence on 
the device characteristics is desirable. 
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